Open system dynamics with non-Markovian quantum jumps 
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We discuss in detail how non-Markovian open system dynamics can be described in terms of 
quantum jumps [J. Piilo et al, Phys. Rev. Lett. 100, 180402 (2008)]. Our results demonstrate that 
it is possible to have a jump description contained in the physical Hilbert space of the reduced 
system. The developed non-Markovian quantum jump (NMQJ) approach is a generalization of the 
Markovian Monte Carlo Wave Function (MCWF) method into the non-Markovian regime. The 
method conserves both the probabilities in the density matrix and the norms of the state vectors 
exactly, and sheds new light on non-Markovian dynamics. The dynamics of the pure state ensemble 
illustrates how local-in-time master equation can describe memory effects and how the current state 
of the system carries information on its earlier state. Our approach solves the problem of negative 
jump probabilities of the Markovian MCWF method in the non-Markovian regime by defining 
the corresponding jump process with positive probability. The results demonstrate that in the 
theoretical description of non-Markovian open systems, there occurs quantum jumps which recreate 
seemingly lost superpositions due to the memory. 

PACS numbers; 03.65.Yz, 42.50.Lc 



I. INTRODUCTION 

The theory of open quantum systems describes the dy- 
namics of a system of interest interacting with its envi- 
ronment The system-environment interaction leads 
to non-unitary reduced system dynamics and the system 
state is described by a density matrix instead of a sin- 
gle state vector used for closed systems. Generally, the 
density matrix evolution is governed by a master equa- 
tion whose unitary part contains the dynamics as given 
by the system Hamiltonian and the non-unitary dissipa- 
tor describes the effects that the environment has on the 
system. 

The presence of the environment leads to decoherence 
which is harmful for practical applications like quantum 
information processing . On the other hand, decoher- 
ence has a role in open fundamental problems of quantum 
physics such as quantum to classical transition [3,] . Often, 
the environment is seen to have unavoidable effects on the 
system dynamics. However, the recently developed abil- 
ity to control quantum systems and the implementation 
of reservoir engineering techniques are revising the role 
of the environment [1, H, Q . This may lead to new ways 
to control the system of interest indirectly via the control 
of the system-reservoir interaction and the properties of 
the environment. 

In memoryless Markovian open systems, the environ- 
ment acts as a sink for the system information. Due to 
the system-reservoir interaction, the system of interest 
loses information on its state into the environment, and 
this lost information does not play any further role in 
the system dynamics. However, if the environment has a 
non-trivial structure, then the seemingly lost information 



can return to the system at a later time leading to non- 
Markovian dynamics with memory. This memory effect 
is the essence of non-Markovian dynamics. 

Non-Markovian systems appear in many branches of 
physics, such as quantum optics [E 0, sohd state 
physics 0, quantum chemistry [lO |. and quantum in- 
formation processing Recently, non-Markovian 
features has also been exploited in the context of 
biomolecules where the environment consists of pro- 
tein solvents However, the elusive nature of non- 
Markovian dynamics makes it often difficult to obtain 
insight into microscopical physical processes governing 
the time evolution. At the same time the complex math- 
ematical structure of the non-Markovian models prevents 
generally to solve the dynamics of the system of interest. 
Hence, new ways to describe non-Markovianity and new 
methods to solve non-Markovian dynamics are highly de- 
sired. 

The density matrix can also be seen as a collection, or 
ensemble, of state vectors. Then, the interaction between 
the system and the reservoir removes the precise informa- 
tion about the specific state vector to describe the system 
state. Instead, the state of the open system is associated 
with an ensemble of state vectors where each state vector 
has a certain (classical) probability of appearance. This 
view has led to the development of Monte Carlo simula- 
tion methods for MarkovianJllJ3JllJIl, IH [11 and 
non-Markovian [H, [2l|, iJlHlilliJopen systems. 
In these methods, the time evolution of each state vec- 
tor in the ensemble contains a stochastic element which 
can be discontinuous (quantum jump) [l^, [3, [H, [H, 
[23, [U, m, m m or continuous (quantum state diffu- 
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One of the most common methods to treat Markovian 
dynamics is the Monte Carlo wave function (MCWF) 
method which exploits quantum jumps ^sj. How- 
ever, a generalization of this Markovian method to non- 
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Markovian regime has turned out to be a challenging 
problem. The central obstacle has been the appear- 
ance of negative quantum jump probabilities due to 
the temporarily negative decay rates of non-Markovian 
dynamics. Earlier approaches to this problem exploit 
auxil iary extensions of the Hilbert space of the sys- 
tem llOl . [20I [2ll . [2^ or exploit the state of the total sys- 
tem I23|. 



We have recently shown that the jump-like unravel- 
ling of non-Markovian master equations is possible within 
the Hilbert space of the system, and hence the auxiliary 
extension of the system Hilbert space is not necessar- 
ily needed [13] ■ The key feature of the developed non- 
Markovian quantum jump (NMQ J) method is the notion 
that, when the decay rates appearing in the master equa- 
tion become negative, the direction of the information 
flow between the system and the reservoir gets reversed. 
During the initial positive decay region, the information 
flows from the system to the environment, while during 
the negative decay the system may regain some of the 
information it lost earlier. In terms of quantum jumps 
this means that the seemingly lost superpositions in the 
ensemble can be restored. This leads to new insight into 
the concept of memory, which is the central ingredient 
of non-Markovian dynamics. We also describe in detail 
the positive and negative factors affecting the numerical 
performance of the method. The ultimate limit for the 
numerical performance is given by the effective ensem- 
ble size Ncs (Sec. IV E[) since the method needs to evolve 
simultaneously A'cff state vectors. 



Our results help to explain why local-in-time master 
equations |26| can indeed describe systems with mem- 
ory and the results also show the presence of some coun- 
terintuitive features of non-Markovian dynamics. In this 
regime, the rate of the process is proportional to the tar- 
get state, instead of the source state, and hence chal- 
lenges the classical view. We show here two different 
proofs of the equivalence between the algorithm and the 
master equation, discuss in detail how the method works, 
and apply it to multi-level atom schemes. Recently, the 
existence of a measurement scheme interpretation of non- 
Markovian dynamics has been actively discussed [13, . 
Our results align along the results of Ref. [28] . We discuss 
this and other insight provided by the NMQJ method in 
detail. 



We have organized the paper in the following way. 
Sec. in] describes briefly the Markovian MCWF method 
and sets the scene for its non-Markovian generalization 
which is presented in Sec. IIIII We then present several 
examples on the use of the NMQJ method in Sec. lIVI and 
discuss the insight provided by the method in Sec.|Vl Fi- 
nally, Sec. IVII concludes the paper. 



II. MARKOVIAN MONTE CARLO WAVE 
FUNCTION METHOD 

Our non-Markovian quantum jump method generalizes 
the MCWF method llSj] into the non-Markovian regime. 
The algorithms and the proof of correspondence with the 
master equation for the two methods are very similar. 
The essential difference is the form of the jump operators 
and jump probabilities. We present first the central in- 
gredients of the Markovian MCWF method and illustrate 
the problems that prevents its use for non-Markovian sys- 
tems. 



A. The algorithm and equivalence with mcister 
equation 

The MCWF method is probably the most commonly 
used Monte Carlo method to treat Markovian open sys- 
tems whose dynamics is g overned by the master equation 
in the Lindblad form [3 (1^ 

3 

Here, p is the density matrix of the reduced system, Hs 
the hermitian system Hamiltonian, Tj is the positive and 
constant decay rate to decay channel j, and Cj are the 
Lindblad (jump) operators describing the effects of the 
environment on the reduced system. 

To unravel the master equation ([T|), MCWF method 
generates an ensemble of stochastic state vector realiza- 
tions whose deterministic and continuous time evolution 
is interrupted by randomly occurring discontinuous quan- 
tum jumps. The average over the ensemble of stochastic 
realizations gives the properties of the reduced system at 
any given moment of time. A generic way to write the 
density matrix in terms of the ensemble is 

pW = E^I^-W)(^"WI' (2) 

a 

where Na{t) is the number of ensemble members in the 
state \'4)a(t)) at time t and N is the total number of state 
vectors in the ensemble (ensemble size). 

The method proceeds in discrete time steps 5t, and 
we consider one step that takes us from time t to t + St. 
During this time step, a given state vector \^ait)) evolves 
either in a deterministic way or performs a randomly oc- 
curring quantum jump. The deterministic evolution is 
given by the non-Hermitian Hamiltonian 

H = Hs~'-^Y.^,C]C,. (3) 
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The essential feature here is the second term on the r.h.s., 
which is constructed from the jump operators that ap- 
pear in the master equation This term reduces, in 
the Markovian case, the occupation probabihty of the 
states which decay. The deterministic time-evolution by 
the Hamiltonian ([3]) leads, for small enough time step St, 
to the state 



\(f>ait + 5t)) 



1 



iH6t 



(4) 



Before the next time step, this state is renormalized and 
the time evolution of \ipa) is 



\\Mt + st))\\' 



(5) 



If, instead of the deterministic evolution, a quantum 
jump to channel j occurs, the state vector changes in a 
discontinuous way 



\cMm\' 



(6) 



The probability for a state vector Itjja} to have a quan- 
tum jump to channel j is directly proportional to the cor- 
responding decay rate Tj, the time step size dt and the 
occupation probability of the decaying state 



(7) 



The choice between the deterministic and jump evolu- 
tions, Eqs. (O and (O respectively, is done by comparing 
a generated random number ^ to the total jump probabil- 
ity Pa . This is the sum over channel specific probabilities 



Pa 



(8) 



and has a direct relation to the norm of \(f)a{t + St)): 

l-pa = mait + 5t))\\\ ^ _ 

By calculating the average evolution cja of |V'q(^)) over 
the deterministic and jump paths one obtains 

\(l)a{t + St)){(j)a{t + St)\ 



aa{t + St) = (l-pa) 



1 -Pq 



{Mt)\C}Q\Mt)) ' 

Here, {I — Pa) is the no-jump probability which weights 
the deterministic evolution and jump probabilities p^, 
weight the corresponding jump paths. 

By inserting Eqs. Q and (O into the Eq. © and rear- 
ranging the terms, one obtains after straightforward cal- 
culation the master equation ([1]) for state vector |?/'Q,(i)). 
Taking a further step by considering the average over the 
whole ensemble, 



a{t + St) 



St), 



(10) 



it is straightforward to see that the master equation ([T]) 
and the MCWF method result given by Eq. pU|) match, 
and the two approaches are indeed equivalent descrip- 
tions of the Markovian open system dynamics. 



B. Why the MCWF does not work for 
non-Markovian systems? 

In Markovian systems, the decay and decoherence pro- 
cesses occur at constant positive rates [c.f. Eq. ([1])]. This 
indicates constant flow of information from the system 
to the environment before the steady state is reached. 
For non-Markovian systems, the decay rates are time- 
dependent and may acquire temporarily negative values 
(to be described in detail in the next Section). During 
the initial period of positive time dependent decay, the 
rate of the information flow changes but the direction of 
the flow remains constant, i.e., from the system to the en- 
vironment. When the decay rate becomes negative, the 
direction of the information flow is reversed and the re- 
duced system, due to the non-Markovian memory, begins 
to recall the information that was lost earlier. 

In the MCWF method, the quantum jump probability 
is directly proportional to the decay rate [c.f. Eq. ([7])] 
which acquires negative values in the non-Markovian 
case. As a consequence of these two facts, a quantum 
jump has negative probability to occur while the deter- 
ministic evolution has larger than 1 probability. There- 
for, it is impossible to make a decision between these two 
alternatives and as a consequence, MCWF method can 
not be used to describe non-Markovian dynamics. 

Earlier attempts to solve this problem exploit usually 
the idea that non-Markovian dynamics can be converted 
to Markovian one by extending the Hilbert space of the 
system [l^ [l^, Hll, [2^ . This may come with a cost for 
computational efficiency and may also prevent obtaining 
insight into non-Markovian dynamics. This also leaves 
open a fundamental question: Is there a corresponding 
jump process in the Hilbert space of the system which 
has a positive probability? 



III. NON-MARKOVIAN QUANTUM JUMPS 

Our starting point is the general local-in-time non- 
Markovian master equation [l], [2^ 

m - l[i/s,p(i)]+E^w^^-w^w^iw 



i^A,(0{p(t),Cj(0Q(t)} 



(11) 



The difference, when comparing to the Markovian master 
equation ([T]) , is that the decay rates A^- (t) depend on time 
and may acquire negative values. In the most general 
case the Lindblad operators Cj{t) may also depend on 
time. 
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A{t) > 



N 



N 




FIG. 1; (Color online) Initially all the A'' ensemble members 
share the same initial state \tpo), i.e., A^o(O) — N. Quantum 
jumps during the positive decay rate (arrows to the right) 
spread the ensemble members to a wider set of different states. 
On the contrary, the non-Markovian quantum jumps during 
the negative decay rate (arrows to the left) transfer ensemble 
members always to states which already exist in the ensemble. 



Special case: Non-Markovian time scale is the 
shortest one 



Before going to the general solution in the next sub- 
section, we first describe the method for the simple case 
in which the non-Markovian time-scale is the fastest one, 
which is most often the case. This allows to introduce 
the NMQJ method in a way that is conceptually rather 
straightforward. With this approximation, the state vec- 
tors do not have time to evolve due to the system Hamil- 
tonian Hs on the time scale of non-Markovian dynamics. 
Consider now a non-Markovian system where the decay 
rates oscillate between positive and negative values be- 
fore reaching a constant Markovian value. For the sake of 
simplicity, we assume here first that all the decay chan- 
nels take negative values simultaneously. At the end of 
the first positive decay period, the initial pure state has 
evolved to a mixed state which can be described in terms 
of the jump paths, unravelled by the MCWF method in 
the positive region, as 



N 



(12) 



Here, lipoit)) is the deterministic evolution from the ini- 
tial state \ipoiO)) without jumps and Itpj) describe the 
ensemble members that have performed one jump to 
channel j, such that \4>j) = Cj\4>o) /\\Cj\^po)\\. In the 
next term, \4'j.k) correspond to members who have per- 
formed first a jump to channel j and then, further- 
more, a second jump to channel k, so that lipj.k) = 
CkCj\ipo) /\\CkCj\ipo)\\. The rest of the terms go in the 
corresponding way. A'o, Nj and Nj^k are the correspond- 
ing numbers of the ensemble members. 

The central question is now how the ensemble is 
evolved so that the result matches the master equation 
PT|) . The sign change of the decay rate indicates the re- 
versal of the information flow between the system and the 
environment, so that for negative decay the system par- 



tially recovers the information that it lost earlier. This 
restoration of lost information is the essence of the non- 
Markovian memory. In other words, the decoherence that 
occurred in the preceding positive decay region, turns to 
re-coherence in the negative decay region, i.e., the earlier 
effects of decoherence get partially cancelled. 

This leads to the idea that non-Markovian quantum 
jumps, taking place in the negative decay region, can- 
cel the effect of the jumps that appeared earlier in the 
positive decay region destroying quantum superpositions. 
Reverse quantum jumps during negative decay are thus 
expected to counteract prior positive decay jumps. This 
means that in the expansion (jl2[) . state \ipj) jumps back 
to the state iV'o), state \ipj,k) jumps to the state {ipj), and 
so on. The direction of the probability flow gets reversed 
for negative decay region as illustrated in Fig. [1] The 
corresponding non-Markovian quantum jump operators 
are 



(13) 



and so on. The probabilities for the jumps to occur are 



No6t\Aj\{Mt)\C}C,\Mt)) 



P,.k 



N. 



(14) 



Equations (fT3|) and (fT4|) demonstrate that the probability 
for reversing a jump for one particular channel is given 
by the portion of ensemble members that have not yet 
jumped in that channel. The numerator gives the to- 
tal jump probability in the ensemble which is distributed 
equally to those ensemble members which can perform 
the jumps. By doing the reversed jump according to 
Eq. (|13p . the discontinuous history of the ensemble mem- 
ber is preserved. This means that when we are reversing 
a jump, we are not erasing the past. 

To prove that the algorithm matches with the master 
equation, we follow very closely the proof of the MCWF 
method [3| ■ The basic idea is to average over the deter- 
ministic and jump paths in order to obtain an equation 
of motion for the reduced density matrix. Evolving the 
ensemble over time step St, gives 

'^it + st) = eoit) + J2^[e,it) + e,^oit)] 



N 



Here, Ooit) is the contribution arising from the deter- 
ministic evolution between times and t. Qjit) is the 
contribution of the ensemble members who jumped ear- 
lier once to channel j and the jump is not cancelled at 
the current point of time. In 0j_>o(t), there has been 
one jump to channel j and which gets cancelled at the 
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current point of time. The rest of the terms arise corre- 
spondingly. It is worth noting that it is not possible to 
cancel something which never happened. Hence there are 
no jumps which can be cancelled from Qoit) part. Taking 
into account for the appropriate weights and keeping in 
mind the jump operators and probabilities from Eqs. ([T3|) 
and (|14p. these terms can be written explicitly 



en = 



\Mt + st)){Mt + st)\ 



1 + no 

\cbj{t + 5t)){ck,{t + 5t)\ 



1 



e,^o(i) = P,^oD,^oH,{t)){^,{t)\Dl^,. (16) 
Here, the time evolved deterministic states are 



\Mt + 6t)) = (i- 
\Mt + st)) - (1- 



iHsSt 
h 

iHsSt 



■m 



(17) 



and their normalization factors are 

no = ^<5t|A„,(i)|(Vo(i)|CiC™|V^o(t)), 

m 



(18) 



All the rest of the terms follow correspondingly. Using 
Eqs. (fT7| and (fT8|) in Eq. (fT6|) and inserting the results 
into Eq. (|15p gives the master equation pT|) . 



In a multi-channel system, positive and negative chan- 
nels may appear simultaneously. The description above 
contains all the negative channels while the positive chan- 
nels evolve according to the MCWF method. Hence, the 
match between the positive channel dynamics with the 
master equation can be proven along the MCWF proof. 
For the sake of simplicity, we leave the detailed descrip- 
tion of simultaneous positive and negative channels to 
the general treatment presented in the next subsection. 



B. General case 

The simplified case presented in the previous section 
nil Al is now generalized. The simple treatment fails in 
a general case, because it assumes that the jump history 
can be unambiguously reconstructed for each state in the 
decomposition p^ . In general, starting from |'(/'o)(V'o|> 
many different combinations of jumps may lead to iden- 
tical contribution \ipa){'4'a\, and all these states should 
be counted together to form Na. 

As in the Markovian case, we write the density matrix 
in the most generic way 



Pit) - E 



N 



(19) 



The positive and negative decay channels are noted with 
j+ and j_, respectively, while the corresponding decay 
rates are Aj^ (t) > and Aj_ (t) < 0. With this notation 
the master equation (jlip can be written as 



m = ^[Hs,pit)] + Y,A,^it)\c,4t)pit)Clit)-^{p{t),cl{t)C,^ 



(20) 



The deterministic time evolution of the state vectors 
\il)a{t)) occurs as before 



|V'«(t)) ^ \M + 5t)) 



(21) 



where the non-normalized state \(j)a{t + 5t)) has been ob- 
tained with the usual non-Hermitian Monte Carlo Hamil- 
tonian. For the sake of convenience, we write this Hamil- 



tonian separating the positive and negative channels 



ih 



J2A,_{t)Cl{t)C,_{t). 



(22) 



The jump probabilities and the jumps for the positive 
channels j+ follow the MCWF prescription, i.e., 

PiHt) A,4t)dt{^^{t)\Clit)C,4t)\i;„{t)), (23) 
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and 



(24) 



term {ipa'{t)\Cj_ {t)Cj_ Moreover, if there are 

no ensemble members in the target state, Na' = 0, then 
the jump probabihty is equal to zero. 



correspondingly. 

For negative channels j- 
cess gets reversed 



The sign of the decay rate Aj{t) can 

when for 



First, 



be understood 
a given chan- 



the direction of the jump pro- 



iiQ-K'W>ir 



(25) 



way. 

0, the process goes as |V') IV'') = 
Later on, when the decay rate be- 
Aj(t) < 0, the direction of this process 
the jump occurs to opposite direction 



In other words, 
takes the form 



the jump operator for negative channels 



i?i-_„,(t) = |^a'(i))(V'o(i)l, 



(26) 



in the following 
nel j, Aj{t) > 

comes negative, 
is reversed and 

Generally, Eq. indicates that the explicit target 
state \i{}a'{t)) of the reverse jump for the source state 
IV'a(t)) is not necessarily unique. This means that the en- 
semble members in the state |V'a(i)) can jump to different 
target states along Eq. (|25p whenever the corresponding 
jump probability is larger than zero. The major factor 
for the computational cost is defined by how many differ- 
ent types of states vectors are created during the positive 
decay region and the need to evolve them simultaneously 
due to their dependence in the negative decay region. 
This point is discussed more in Section IV El 

The proof of our NMQ J method follows again the same 
lines of the Markovian MCWF method ^IS*! and given in 
the previous Section. By weighting the deterministic and 
(27)jump paths over the time step St with the appropriate 
probabilities we obtain the master equation pT|) . Calcu- 
Note that the probability of the non-Markovian jump lating the average a of the evolution of the ensemble (fTQ]) 
is given by the target state {ipa') of the jump along the over St gives 



state of the jump is lipait)) = 
_{t)\ipct'it))\\- The source and target 



where the source 

Q_(t)|Vo' (<))/! IQ 

state of the jump swap their role when the decay rate 
becomes negative. 

This transition for a given state vector \ipa) hi the en- 
semble (|19|) occurs with probability 



Na,{t) 
N^t) 



A,_ {t)\St{^^,{t)\d_ {t)C,_ (i)IV'a'(O)- 



CT{t + St) J2 



Nqjt) 

N 



3+ 



\Mt + st)){Mt + st)\ 
\\\Mt + st)W 



„ . C,^(t)\lPJt))(lpait)\C] (i) _ . - , 

^""""^^^^^JWW^^ ^ e-.«'W^i--..'WIV'.(t))(V'«(t)|i^i-J„,(i) 



(28) 



Here, the summations a and a' run over the ensemble 
[c.f. Eq. (jl9[) ]. the summation over j+ and j_ cover the 
positive and negative channels, respectively. The first 
term on the r.h.s., in the summation over a, is the prod- 
uct of the no-jump probability and the deterministic evo- 
lution of the state vector, the second and third terms de- 
scribe the positive and negative channel jumps, respec- 
tively, with the corresponding probabilities. 

The details of the proof are presented in Appendix A 
and we describe here briefly the main features. Like in 
the Markovian MCWF case, the deterministic evolution 
gives the commutator and the anticommutator parts of 
the master equation. Moreover, the jump part of pos- 
itive channels goes along MCWF giving the remaining 
"sandwich" term for positive channels After making 



the series expansion of the denominator of the determin- 
istic part and keeping the terms to the first order in St, 
we are left with the norm change term due to negative 
channels times the deterministic evolution, jump proba- 
bility for negative channels times the deterministic evolu- 
tion, and the jump term for negative channels. As shown 
in Appendix A, the first and last of these three cancel 
each other and the second one gives the " sandwich" term 
of the master equation ([^0]) for negative channels. This 
completes the proof. 



IV. EXAMPLES 

In order to demonstrate the appUcabihty of the NMQJ 
method we give now concrete examples. These examples 
also show how the method works at the level of single re- 
alizations in the ensemble. Our physical system of choice 
is an atom interacting with a Lorentzian structured reser- 
voir, e.g., an atom interacting with a single mode of a 
leaky cavity. 

The first example is a two-level atom interacting off- 
resonantly with the cavity field, also known as detuned 
Jaynes-Cummings model [c.f. Fig. [2] (a)]. We use this 
simple system to give a detailed walk-through description 
on how the NMQJ method is implemented in practice. 
The other examples deal with a three-level atom, another 
archtype of atomic systems, which holds two independent 
decay channels and also three different level geometries: 
A, V, and ladder-systems [c.f. Fig. [2] (b)-(d)]. For these 
cases we see how having simultaneously both a negative 
and a positive channel results in rich dynamics. 

The structure of the effective ensemble, i.e., states 
and the way in which they connect by different jump 
channels, is shown in Fig.[3]for each example case, respec- 
tively. This illustrates, how physically identical states 
can be reached by different combinations of jumps in the 
V and ladder-systems. 

From the NMQJ method's point of view the details of 
the actual physical system and the variety of approxi- 
mations during the derivations are irrelevant as long as 
the master equation is in the desired general form, given 
by Eq. pT|) . Moreover, just to highlight this feature, we 
illustrate explicitly how the NMQJ method follows the 
formal mathematical solution of the given master equa- 
tion as far as the solution is physically consistent, i.e., the 
density matrix remains positive. If the solution fails to 
be positive at some point, it obviously means that some 
of the approximations made while deriving the master 
equation of the reduced system are not valid. 



A. Derivation of the non-Markovian local-in-time 
Master equation 

To give an idea of how non-Markovian local-in-time 
master equations can be derived microscopically and of 
the explicit form of the time-dependent decay rates, we 
give a brief sketch of the derivation for the example sys- 
tem in hand. 

The system Hamiltonian of a multi-level atom is 

i 

Similarly, the self-Hamiltonian for the electromagnetic 
field constituting the environment is 

^^onv = ^ hvkalak- (30) 

k 
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(a) (b) (c) (d) 




FIG. 2: Example cases introducing the level notations, (a) 
Jaynes-Cummings model, (b) A-system, (c) V-system, and 
(d) ladder-system. Only transitions expressed by arrows con- 
tribute since they reside close to the resonance frequency of 
the cavity. 



(a) (b) (c) (d) 




FIG. 3: Effective ensembles for the example cases. (a) 
Jaynes-Cummings model, (b) A-system, (c) V-system, and 
(d) ladder-system. The number in the arrow indicates the 
jump channel. The expressions for states \'4!a) are given in 
the text. 

The dipole interaction between the system and its envi- 
ronment is described by an interaction Hamiltonian 

i/int = -D • E, (31) 

where D = gr is the dipole moment operator and E the 
quantized electromagnetic field. Within the second or- 
der time-convolutionless (TCL) approach [ij and after 
performing the secular approximation, the jump chan- 
nels are categorized by atomic transition frequencies, or 
Bohr frequencies, lo, such that the Lindblad operators 
are 

OJj —UJi—U! 

where dy = {i\{—'D)\j)/D is the dimensionless value of 
the matrix element of the dipole moment operator D 
(dimensional unit D). It is convenient to pass to the 
continuum limit of environmental modes Vk such that 
J2k lo^feP ~^ /di^ J(j^)- Here, at describes the coupling 
strength between the system and the reservoir mode 
j/fe, and J(i^) is the spectral density of electromagnetic 
modes 1]. Considering only the zero temperature envi- 
ronment, where all the modes are initially empty, each 
decay channel is related to a time-dependent decay rate 

poo 

A^{t) = 2 ds di^J{i^)cos[{iy-uj)s]. (33) 
Jo ^0 

The interaction with the reservoir introduces a renor- 
malization of the system Hamiltonian Hs by a Hermitian 
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term, i.e., the Lamb shift Hamiltonian 

HLs{t)^nY,^Ut)ClC^, (34) 

where the time-dependent rate factor is 

X^{t)=l ds di^J{i^)sm[{iy-Lj)s]. (35) 
Jo Jo 

We label the different Bohr frequencies by {uj-'}, where 
j = 1, . . .. Correspondingly, the jump operators are 
Cj = Ci^j and the decay rates are Aj (t) = A^j (i) and 
Xj(t) = X^^j(t). Then, the time- local master equation in 
the interaction picture is in the form of Eq. (Ilip , where 
system Hamiltonian Hs has been replaced by HLs{t). 

The spectral density of the electromagnetic field inside 
an imperfect cavity is well approximated by a Lorentzian 
distribution 



Jlo 



a 



2-K {v 



(r/2)2 



(36) 



where is a coupling constant, Wcav is the resonance 
frequency of the cavity and T characterizes the width of 
the distribution. The essential parameter in this case is 



the detuning 5j 



U!^ of the Bohr frequency with 



respect to the cavity resonance frequency. 

Since the cavity supports only modes residing close 
to its resonance frequency cjcav, only transitions whose 
Bohr frequencies are close to this value contribute to the 
dynamics. This justifies the description of the atom's 
Hilbert space consisting effectively of only two or three 
levels, which we now study. 



B. Units and parameters 

In the examples, the time scale is set by the inverse of 
the spectral distribution width F^^. The resonance fre- 
quency is assumed to be large Wcav ^ T- The Markovian 
time scale is then tm ~ 10 F"^ [c.f. convergence of the 
decay rates to steady Markovian values in, e.g., FigdJ^a)]. 
In the Jaynes-Cummings model the coupling constant is 
set to = 5 and in the three-level systems it is = 2. 
The dipole moment matrix elements are always assumed 
to be dij = 1 for all pairs of states i ^ j. In the numeri- 
cal simulations the time step size is 6t = O.OIF"^ and the 
size of the ensemble is = 10^. The notation of atomic 
levels is the same as in Fig. [21 



C. Results 

For the sake of comparison, we solve the master equa- 
tion in two different ways. First, we solve the density ma- 
trix by using the NMQJ method. Second, we calculate 
the formal analytical solutions of the equations of motion 
of the individual density matrix components (expressions 
are given in Appendix B). The results are then compared 
in order to verify the functionality of our method. 



1. Two-level atom: detuned Jaynes-Cummings model 

The two-level case involves only one Lindblad operator 
Ci — (T„ = which is the usual lowering operator 

from the excited to the ground state. We choose the de- 
tuning (5i = 5 F and the Fig. |4] (a) shows the oscillatory 
behavior of the corresponding decay rate Ai{t) . The ini- 
tial state is a pure state p(0) = \ijjo{0)){ipo{0)\, meaning 
that all the ensemble members are initially in the same 
state |i/'o(0)). In our example |'0o(O)) = (3|a)+2|6))/\/T3. 

For the given single jump operator and an initial state 
including a finite excited state component, there will be 
only two kinds of states contributing to the master equa- 
tion solution. This is because according to the unraveling 
in Eq. ^ the global phase factors of the single ensemble 
members do not affect the density matrix representation. 
The two non-equivalent states are now the evolved ini- 
tial state vector \ipo{t)) and the ground state = \b), 
which can be reached from \tpo{t)) by operating with the 
Lindblad operator. Correspondingly, there are two dis- 
crete variables NQ{t) and Ni{t) counting the number of 
ensemble members on each of these two states. Initially 
iVo(O) = N and iVi(O) 0. 

For a certain initial time interval, the decay rate Ai 
is positive (see Fig. During this period the ensem- 
ble evolves according to the standard MCWF descrip- 
tion. The deterministic evolution \ipait)) \ipa{t + 5t)) 
is given by Eq. ([2^ with Hamiltonian H = H^s — 
f Ai(i)C|Ci = h[Xi{t) - |Ai(t)]|a)(a|. The deter- 
ministic evolution is interrupted by quantum jumps 
\ipo{t)) \^i{t)) occurring with a probability Pgit) = 

A,it)St{Mt)\ClCi\Mt)) = A,it)St\{a\Mt))\'' given 
by the Eq. (I23p . In our notation this means that when 
quantum jump occurs, the occupation numbers are up- 
dated as {No{t),Ni{t)} -> {No{t) - l,Ni{t) + 1}. Once 
an ensemble member has jumped to the state it 
can not experience any other quantum jumps during 
this period, since the corresponding jump probability is 
Pi cx |(a|Vi)P =0. 

After the first positive period the decay rate becomes 
negative. The deterministic evolution is still driven by 
the same Hamiltonian as previously. However, now 
those ensemble members which had previously jumped 
to the ground state \ipi) are able to make a reverse non- 
Markovian quantum jump |'0o(i)) IV'i) going back to 
the dcterministically evolved initial state. The probabil- 
ity of this jump is given by Eq. (|27p and is 



^\A,{t)mMt)\clc,\Mt)) 

= ^\A,{t)\St\{a\Mm'- (37) 

Accordingly, the occupation numbers are updated after 
each reverse jump such that {No{t), Ni{t)} {No{t) -\- 
1, Ni{t) — 1}. The ensemble members in the state \ipo) are 
not able to perform quantum jumps during this period, 
since in the ensemble there are no states \ipa) for which 
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FIG. 4: (Color online) Dynamics of the Jaynes-Cummings 
model. Initial state is \ipo{0)) = (3ja) + 2\b))/Vl3. (a) Decay 
rate Ai(t). (b) Populations paa (initially higher line) and ptb 
(initially lower line), (c) Absolute value of the coherence pab- 



\Mt))^cM/\\cM\\. 

The successive periods of positive and negative decay 
rate are treated in a similar way. In the Fig. [3] we show 
how the ensemble average of single realizations generated 
by the NMQ J method gives the exact solution of the mas- 
ter equation. In the corresponding Markovian case with a 
constant decay rate AMarkov = limt^oo A(t), the solution 
would be a simple exponential decay towards the ground 
state accompanied by exponential decoherence. The non- 
Markovian time-dependent decay rate leads to a slower or 
faster decay compared to the Markovian exponential one. 
Furthermore, since the decay rate takes negative values, 
the decay process can be partially reversed. This leads 
to a regain of excited state probability and re-coherence. 

In Fig. [5] we give an example of a single realization 
experiencing both a quantum jump to the ground state 
\ipo) — > lipi) during the positive decay rate and a reverse 
non-Markovian quantum jump back to the initial state 
\ipo) <— \ipi), during the negative decay rate. The essence 
of this illustration is that after these two jumps the state 
is (up to an irrelevant global phase factor) precisely the 
same as if the evolution would have been purely deter- 
ministic. However, when evaluating the time evolution 
with the ensemble average, the total contribution of this 
realization is different from the contribution given by a 
realization with no jumps. 
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FIG. 5: (Color online) Example of the dynamics of a single 
realization in the Jaynes-Cummings model (parameters as in 
Fig. m . In this realization the deterministic evolution of the 
initial state \ipo) (higher line) was followed until t « 0.4 F"^, 
when a Markovian quantum jump (arrow down) brought the 
state to the ground state = \b) (lower line). Then, de- 
terministic evolution of the ground state was followed until a 
non-Markovian quantum jump (arrow up) recreated the |?/;o) 
state aX t 0.8 F~^, whereafter the evolution was determin- 
istic. 



2. Three-level atom: K-system 

In a A-system there are two jump channels with Lind- 
blad operators Ci = \h){a\ and C2 = \c){a\. In our 
example we choose the corresponding detunings to be 
5i — — 3r and 82 = 5T. With these values the two 
decay rates have at certain time intervals opposite signs 
[c.f. Fig.[6](a)]. We now look at the initial state \iPq{0)) — 
(4|a) + 2|6) + |c))/V2T. 

Starting with such an initial state the ensemble consists 
of effectively three different states: \4>o{t)), \^pi) = 
and IV'2) = |c). There are now two competing processes 
affecting the time evolution of the initial state. Initially 
both decay rates are positive, but at i w 0.5 F^^, channel 
2 becomes negative. This means that after this moment, 
on the one hand, there are still quantum jumps through 
channel 1 away from the initial state \ijjo) li/ji), but on 
the other hand, channel 2 repumps the ensemble mem- 
bers back to the initial state by non-Markovian quantum 
jumps iV'o) ^ IV'2)- At i « 1.2 F~^ both decay rates 
change their signs and so on all the way until t « 2.5 F~^. 
Fig. [6] illustrates that when the decay rates are counter- 
acting each other, plateaus in the evolution of the density 
matrix elements can be observed. 



3. Three-level atom: V-system 

In the case of a V-system, the two jump channels are 
Ci — \c){a\ and C2 = \c){b\. We choose the detunings as 
earlier: Si — — SF and $2 = 5F. We consider the initial 
state |-0o(O)) = (|a) |6) |c))/V3. In this case, the 
ensemble consists of effectively only two different states. 



10 




FIG. 6: (Color online) Dynamics of a A-system with an ini- 
tial state \ipo{0)) = (4|a) + 2\b) + \c))/V21. (a) Decay rates 
have momentarily opposite signs, (b) A plateau emerges to 
the excited state population paa (initially highest line) as the 
decay channels counteract each other. Other populations ptb 
(initially middle line) and pcc (initially lowest line) behave 
according to the two separate decay channels, (c) Also coher- 
ences pab, Pac, and pcb (initially highest, middle, and lowest 
line, respectively) have plateaus. 



FIG. 7; (Color online) Dynamics of a V-system with an initial 
state IV'o(O)) = (|a) + \b) + \c))/V3. (a) Decay rates, (b) 
Populations of the upper states paa and pbb (lowest and middle 
line, respectively) decay according to single decay channels 
and the ground state population pcc (highest line) increases 
correspondingly, (c) The upper state coherence pab (lowest 
line) has a plateau as the decay channels are counteracting 
each other, while coherences involving the ground state pac 
(middle line) and pbc (highest line) are related to individual 
decay channels. 



since both Lindblad operators act as \ipo{t)) — > jf/'i) = 

The dynamics in Fig. [7] shows how the upper state 
probabilities decay according to the individual decay 
channels. Since only lipo) carries coherences, there is a 
plateau in the upper state coherences, as it is affected 
simultaneously by decoherence and recoherence. 



4-. Three-level atom: Ladder-system 

The ladder-system induces the most complicated dy- 
namics of the three three-level atomic schemes consid- 
ered here. The Lindblad operators form a short cascade, 
Ci = \b){a\ and C2 — |c)(6|, so that the target state of the 
upper channel can still decay further by another quan- 
tum jump. There are three different possible quantum 
jump processes: |V'o(i)) \tpi) = |^) through channel 1, 
\ipo{t)) — > 1^2) = |c) through channel 2, and \tpi) \1lJ2) 
through channel 2. Therefore, the effective ensemble con- 
sist of three state vectors. 

During the negative period of channel 2, there are 



now interestingly two possible target states for a non- 
Markovian quantum jump from state \ip2) correspond- 
ing to processes |V'o(*)) ^ IV'2) and ^ \^2)- The 
example dynamics in Fig. [5] shows how the initial state 
l^o(O)) = (4|a) + 2\b) + |c))/V2T evolves. It is evident, 
that eventually the state decays towards |i/'2), but due 
to complicated connections between the states and the 
changing signs of the decay rates, the dynamics is more 
rich than in the other cases. 

Starting from an initial state \ipoiQ)) — \a}, our other 
example of ladder-system dynamics shows that the den- 
sity matrix loses its positivity at t w 1.0 F"^ (c.f. Fig. [5]), 
which indicates that the approximations in the derivation 
of the master equation do not hold for this level geome- 
try. The NMQ J solution follows the formal mathematical 
solution as long as it remains positive and the method is 
able to identify the point where the time evolution be- 
comes unphysical. The failure of the positivity occurs, 
when channel 2 is still negative while all the ensemble 
members in state IV'2) have already had a non-Markovian 
quantum jump to states \4>o{t)) and IV'i)- This happens. 
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FIG. 8: (Color online) Dynamics of a ladder-system with an 
initial state \ipo{0)) = {4\a) + 2\b) + \c))/V21. (a) Decay 
rates, (b) The ladder structure is clearly visible as the decay- 
ing population of the highest excited state paa appears first 
as increase of the middle state population ptb, and eventually 
everything ends up to the ground state pcc (initially highest, 
middle and lowest line, respectively), (c) Only the initial state 
IV'o) contributes to the coherences pab, pac, and pbc (initially 
highest, middle, and lowest line, respectively), and therefore 
their dynamics is as simple as in the V-system. 



because the probability for such a non-Markovian quan- 
tum jump is Pi^^ oc Na{t)/N2it), where N2{t) 0. 
This property has some interesting imphcations in the 
search for a positivity conditions for non-Markovian sys- 
tems □01 . 



V. DISCUSSION 

A. On non-Markovian quantum jump operators 
and probabilities 

To circumvent the problem of the negative probabil- 
ities of the Markovian MCWF method, one is tempted 
to consider negative probabilities as positive ones for in- 
verted jumps, i.e., to switch the role of the initial and final 



states of a given Lindblad operator by setting Cj 



'0 — J- 

However, this does not lead to the correct ensemble for 
the non-Markovian dynamics. The essence of the nega- 
tivity of the decay rate is the reversal of the decoherence 
process, i.e., re-coherence, and partial cancellation of the 
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FIG. 9: (Color online) Dynamics of a ladder-system starting 
from an initial state |^o(0)) = \a); other parameters are as 
in Fig. [S] Populations of the middle state pbb (higher line) 
and the ground state pcc (lower line) are shown. The NMQJ 
solution follows the formal analytical solution of the master 
equation faithfully until the state loses its positivity at t « 
1.0 as the ground state populations tends negative (thin 
black line marks the zero value). 



decoherence which occured in the past. If one uses in 
the non-Markovian region with negative decay rates the 
substitution Cj Cj , this only replaces one decoherent 
process with another one. 

Let us illustrate this with the simple example we con- 
sidered in Sec. IIV C II Writing the equations of motion 
explicitly for the density matrix elements of a two-level 
system, gives for the positive decay rate region. 



Paa = -|A|pa 

Pbb = IAIpqq, 



Pab 



1, A , 

-l^\'^\Pab, 



(38) 



and for the negative decay region 



A|pa 
-|A|paa, 

Pab = l^\'^\Pab- 



Paa — I^IPaa; 
Pbb 



(39) 



Here, a denotes the excited and b the ground state of 
the two-level atom. The first line of Eq. ([55]) shows that 
during the negative decay the excited state probability 
increases and that this increase is directly proportional 
to the probability which the excited state already has. 
This is a counterintuitive feature since it means that the 
total rate of the process is proportional to the target 
state, and not to the source state as in the positive decay 
region [c.f. Eq. The last fine of Eq. ^ shows that 

the coherences increase during the negative decay. 

If one attempts to remedy the negative probability of 
the jump given by the Markovian method by changing 
the sign of the decay rate and substituting Cj C|, or 
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cr_ — > cr+, this gives the equations of motion 

Paa = |A|pf,f,, 
Phh = -|A|p6b, 

pah - -^|A|Pa6. (40) 

It is easy to see that these equations are not the correct 
equations of motion p9p . In particular, the proportion- 
ahty of the state populations for paa and go wrong, 
and the coherences decrease while the correct equations 
([55)1 show that they must increase. 

Generally speaking, a simple sign change of the de- 
cay rate from positive to negative in the non-Markovian 
master equation (|20p may seem a 'priori as a rather triv- 
ial problem to solve. However, as the simple example 
above illustrates, the sign change actually leads to a very 
complicated problem. The main source of the complica- 
tion is that the non-Markovian jump operators, given by 
Eq. ([M)) . do not appear explicitly in the master equation 
to be solved, whereas in the Markovian case one can pick 
the jump operators directly from the dissipator of the 
master equation. 

It is also interesting to note that we can interpret the 
jump probability (|27p in the following way. The numer- 
ator 7V„/|Aj_(i)|(5i(V'a'(0|C]_ WCj-WIVJa'W) gives the 
cumulative non-Markovian quantum jump probability in 
the whole ensemble. This is then divided to those iV^ 
ensemble members I'j/'q) who perform the jumps. 

B. Why local-in-time master equation can describe 
non-Markovian dynamics with memory? 

Two common ways to describe non-Markovian open 
system dynamics are the memory kernel master equa- 
tions and the local-in-time master equations with time 
dependent decay rates [l|. The former consists of an 
integro-differential master equation where the change of 
the system state at a given moment of time is given by 
the integral over the past evolution according to a given 
memory kernel. The local-in-time master equations in 
turn are based on the microscopic system-resevoir in- 
teraction modeling leading to a differential equation of 
motion for the density matrix of the system, which is 
local-in-time. The description of non-Markovian dynam- 
ics without the use of a memory kernel, as done with the 
local-in-time master equations, may seem at first sight 
counterintuitive. Our NMQJ method sheds new light on 
this issue and shows explicitly how and where the mem- 
ory appears in local- in-time master equations. 

Suppose now that we have a density matrix of the sys- 
tem p(i) and the corresponding ensemble of state vectors 
during the initial positive decay region. At each time step 
a certain small fraction of the state vectors may jump 
to decay channel m according to the Markovian MCWF 
scheme: -> IV-) = C™| |C™|V'') 1 1- K is impor- 
tant to note that ji/;') contains the information what the 



state IV') was before the jump |i/'') IV') took place, and 
that the whole ensemble still includes both types of state 
vectors 1-0') and |V'). Then the system enters into the neg- 
ative decay rate region. Here, as described in the previ- 
ous two subsections, the jumps go into opposite direction 
from to IV''): and the probability of this jump is given 
by the target state |'0')- In other words, the very state 
vector that contains information on the past state of jV') 
defines both the target state of the non-Markovian jump 
and the probability for this jump to occur. In this way 
the past affects the current evolution of the system [3l|. 

It is difficult to see from the density matrix description 
where the memory of the earlier state of the system is. 
However, according to the description above, when we 
look at the density matrix as an ensemble of state vec- 
tors and study the dynamics of the state vectors in terms 
of the jumps, we see explicitly how the ensemble mem- 
bers carry memory of other ensemble members. This 
memory comes into play when the decay rate becomes 
negative. It is also important to note that if the number 
of ensemble members in the target state of the reverse 
jump becomes equal to zero, Na' — 0, then the system 
has lost its memory, and consequently the reverse jump 
probability vanishes since it is directly proportional to 
N^, [c.f. Eq. 



C. Is continuous measurement of environment 
allowed for non-Markovian systems? 

For Markovian open quantum systems, single Monte 
Carlo realizations have a measurement interpreta- 
tion p^ . The environment is thought to be monitored 
in a continuous way, and the corresponding reduced sys- 
tem evolution, conditioned on the measurement outcome, 
constitutes a single pure state trajectory of the ensemble. 
The existence of a measurement scheme interpretation 
for non-Markovian trajectories has been recently under 
active debate. Diosi claims that, at least in principle, 
certain types of QSD trajectories can be interpreted as 
true pure state single system trajectories [l^]. His idea 
is based on the assumption of availability of an infinite 
set of entangled von Neumann detectors. Wiseman and 
Gambetta question Diosi's claims and the existence of 
true pure state trajectories with the measurement scheme 
interpretation. Their argument is based on the notion 
that in Diosi's scheme one should actually measure also 
those von Neumaim apparatuses which are yet to interact 
with the system 2S]. Due to the entanglement between 
the von Neumann apparatuses, the measurement induces 
noise turning the true pure state trajectories into mixed 
ones. 

Though both works mentioned above deal with diffu- 
sion descriptions, it is interesting to note how our jump 
scheme fits into the discussion. In the NMQJ method, 
the memory of one ensemble member is carried by other 
ensemble members. When a reverse non-Markovian jump 
for a given ensemble member occurs, this member returns 
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to the state which it would have at this point of time, if 
the prior positive decay jump had not occurred. In the 
simple two-level atom example, the superposition which 
was lost earlier gets restored by the non-Markovian jump, 
and the information on the earlier state of the system re- 
turns from the environment to the system. The crucial 
point is that the information lost by the system to the 
environment in the initial positive decay region has to be 
still available to the system when the decay rate turns 
later on negative. If we measure the environment in a 
continuous way, we are extracting information from the 
environment - and indirectly on the system state. If this 
measurement is destructive, then the information is not 
available to the system anymore and the non-Markovian 
dynamics gets distorted. In the case of a two-level atom, 
the measurement of the photon in the environment de- 
stroys the photon, and the two-level atom can not get 
re-excited during the negative decay region. 

In addition, in the two-level atom example, the oscilla- 
tions in the excited state probability arise due to vir- 
tual exchanges of excitations between the system and 
the reservoir P, [2l], [2^. Virtual processes can not be 
directly measured while they still affect the system dy- 
namics. This fits to the insight that the NMQJ gives 
though in terms of virtual processes there is a subtle dif- 
ference: instead of virtual exchange of photons between 
the two-level atom and the reservoir, we rather describe 
the oscillations in the excited state amplitude of the atom 
as destruction and restoration of the quantum superposi- 
tion. This difference between the two descriptions arises 
because an absorption of the photon by the atom means 
a jump from the ground state to the excited state. This 
process, by definition, can not increase the coherences 
which is a key feature of non-Markovian systems in the 
negative decay region, as discussed in detail Sec. IV Al 

If single realizations can not be measured, is there some 
other physical meaning that they have? In our formalism 
the probability to be in a given state at a given moment 
of time is the sum of all the paths leading to this state, 
see Fig. [TOl In this sense the state vector evolutions can 
have an interpretation as possible paths that the system 
may take from its initial to final state. However, com- 
bining with the lack of measurement scheme, this means 
that we are not allowed to measure which path the sys- 
tem has taken while all possible paths contribute to the 
system state. If we try to extract information on the fol- 
lowed path by means of measurements, we disturb the 
non-Markovian memory. The rigorous connection to the 
Hilbert space path integral formalism will be studied in 
the future. 



D. Basic comparison to other jump descriptions 

Earlier approaches to treat non-Markovian dynamics 
with quantum jumps use auxiliary states and exploit the 
idea of Markovian embedding of non-Markovian dynam- 
ics in the extended Hilbert space [H [13, [2l|, [Hi . Other 




► 

TIME 

FIG. 10: (Color online) (a) Sketch of a time-dependent decay 
rate with periods of positive and negative values (arbitrary 
units) . (b) Examples of single realizations encountered in the 
ensemble. The system is assumed to be such that there is only 
one decay channel and two physically different states: the 
initial state \tpo{t)) and the target state of a quantum jump 
IV'i(i)) (deterministic evolution is given by thin horizontal 
lines). The state of an ensemble member at the given time 
is indicated by the thick line. Quantum jumps from \tpo) to 
IV'i) (arrows down) occur at random times during the positive 
decay rate, while non-Markovian quantum jumps from 
back to It/jq) (arrows up) occur during the negative decay rate. 
The total probability to have state {ipo} at the end of the 
shown evolution period is the sum of the paths 1 and 4 while 
the probability to have state \tpi) is the sum of the paths 2, 
3, and 5. 



jumplike unravelings use as an aid the state of the total 
system and hidden variables [22] or take the measure- 
ment theory perspective [32] . Our results show that it is 
possible to have jump-like unraveling of non-Markovian 
dynamics of the reduced system without extending the 
system Hilbert space or considering in detail the total sys- 
tem dynamics and hidden variables. It is worthwhile to 
see if the differences between our method and those devel- 
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oped earlier reveal interesting aspects of non-Markovian 
dynamics. For this purpose, we compare our method 
to pseudomode (PM) [lol, doubled ^| and triple ^ 
Hilbert space methods (DHS and THS respectively) , and 
to the (quantum trajectory method based on hidden vari- 
ables 

The PM method describes the properties of the envi- 
ronment in terms of the auxiliary pseudomode(s) with 
whom the system of interest interacts [l^. The pseudo- 
mode is then coupled to the Markovian reservoir while 
the system of interest interacts only, in a coherent way, 
with the pseudomode. The Markovian pseudomode mas- 
ter equation can be unravelled with the MCWF or some 
other Markovian method. Once this is done, the dy- 
namics of the system of interest is obtained by tracing 
out the pseudomode. This leads necessarily to mixed 
state trajectories for the system of interest while in our 
NMQJ method the time evolution of the ensemble mem- 
bers consists of pure states living in the Hilbert space 
of the system. In addition, the PM method relies on 
some assumptions on the form of the environment spec- 
tral density so that the pseudomode structure can be 
calculated, and it also exploits the solution of the total 
system dynamics. Our NMQJ method differs from the 
PM method in both of these issues and has been used 
to simulate two-level atom in photonic band gap in the 
absence of driving between the two states [l^l (the driven 
case is more challenging, see the next subsection). 

On the other hand, the pseudomodes are by construc- 
tion directly related to the properties of the environment. 
As a matter of fact, it is possible to show by exploiting 
the insight given by the NMQJ method, that the pseu- 
domodes can be interpreted as an effective description of 
the memory of the environment of the open system [ssj . 
This is based on the notion that periods of negativity of 
the decay rate of local-in-time master equation coincide 
with those periods of time during which the pseudomode 
feeds coherently the system. 

The doubled Hilbert space (DHS) method uses two 
copies of the state vector to create a single realization in 
the ensemble The time evolution of the two copies is 
identical in the positive decay region. When the jumps 
with the Lindblad operators occur during the negative 
decay, one of the two copies gets multiplied by —1. This 
produces a negative contribution to the ensemble aver- 
age. The probability in the ensemble is conserved be- 
cause the norm of the deterministically evolving state 
vectors increases to values larger than one. From the 
statistics point of view, this means that the number of 
jumps during negative decay has to match the increase of 
norm in the deterministic evolution, and the probability 
is conserved on average. The consequence is an addi- 
tional source of statistical noise. In the NMQJ method 
each state vector is normalized to one at each time step 
and the probability is conserved exactly. This gives a 
better statistical performance over the DHS method. In 
addition, the NMQJ avoids the numerical burden which 
is present in the DHS method due to the doubling of the 



Hilbert space size 34 1. 

An interesting improvement to the DHS method is pro- 
vided by the triple Hilbert space (THS) method [23| . 
This method shows that the Markovian embedding of 
non-Markovian dynamics can be done with only three 
auxiliary discrete states. The original system dynam- 
ics is then contained in the coherences of the extended 
space state vectors. The method avoids the additional 
statistical noise term of the DHS method. However, the 
THS method uses a 4 times larger number of decay chan- 
nels and a 3 times larger Hilbert space than the NMQJ 
method. Moreover, since the dynamics of the original 
system is contained in the coherences of the extended 
space, unphysical situations such as violations of posi- 
tivity of the density matrix during the time evolution 
may occur and pass unnoticed. In contrast, the NMQJ 
method, by construction, always keeps the dynamics pos- 
itive since it is not possible to a have negative integer 
number of state vectors in the ensemble. It is also worth 
mentioning that in the THS method the auxiliary quan- 
tum jump channels open when the decay rate becomes 
negative. This means that during negative decay inter- 
val, the probability flows out of the Hilbert space of the 
original system whereas in the NMQJ method, the direc- 
tion of the probability flow within the Hilbert space of 
the system gets reversed at this point. 

From the fundamental quantum physics point of view, 
it is also interesting to discuss the jumplike unraveling of 
non-Markovian dynamics which is based on hidden vari- 
ables [l^l . The basic idea of the method is to obtain the 
system trajectories from the guiding state describing the 
state of the total system. This is then used to obtain the 
stochastic evolution of the so called property state which 
includes information on the value of the environmental 
hidden variable and the corresponding properties of the 
reduced system. Our result seems to indicate that it is 
possible to describe non-Markovian dynamics with quan- 
tum jumps without the use of hidden variables. However, 
since the hidden variable approach allows jumps towards 
ground and excited states in the two-level atom case, it 
would be very interesting to compare in detail the time 
evolution of the ensemble members in both of the meth- 
ods, and to see if there exists any connections between 
the two. 



E. Numerical and technical aspects 

Since in the NMQJ method the realizations depend 
on each other due to memory effects [c.f. Eq. (|77|) ]. it 
seems at first sight that all the N ensemble members 
have to be evolved simultaneously. However, according 
to Eq. the ensemble consists of several copies of 

each \ipa{t)). Obviously, there is no need to have on a 
computer several copies of the same state vector. It is 
sufficient to have one copy and the corresponding integer 
number Na- Any number N of the realizations of the pro- 
cess can be done by making N^s ^ N state vector evo- 
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lutions where N^s is equal to the number of terms in the 
summation N = i^-^- When the reahza- 

tions of the process are generated on a computer, a jump 
means changing the integer numbers Na(t) accordingly 
in Eq. ^9]) . A saving in CPU time is achieved since it is 
not necessary at each point of time to evolve N state vec- 
tors, instead, it is enough to decide N times if the jumps 
occurred or not. This means that the NMQJ method 
has a built-in optimization which can be exploited to im- 
prove the efficiency of the method. For the two-level 
atom example described above, the effective ensemble 
size iVeff = 2 while N = 10^. However, these state 
vectors need to be evolved simultaneously since there is 
a dependence between the state vectors [c.f. Eqs. ([25]) - 

m]- 

We can summarize the key factors for the numerical 
performance of the NMQJ method as follows: (i) no 
Hilbert space extensions are needed, (ii) the identification 
of the negative rate process as reverse jumps which keeps 
iVcff constant during the negative decay region and al- 
lows to technical optimization of the simulations (iii) the 
computational cost increases when the number of terms 
in the summation (fTO|) increases. The first two points 
allow to improve the efficiency while the third point is 
expected to set the ultimate limit for the required com- 
putational resources. In addition of this resource limit, 
there exists also non-Markovian systems for which it is 
very challenging to derive local-in-time master equations 
of the form pT|) . An example of this type of the system 
is a driven two-level atom in a photonic band gap ma- 
terial. To the best of our knowledge, there does not yet 
exist local-in-time master equations of the form (lll|) for 
this system. On the other hand, it is possible to simulate 
this system already, e.g., with the method developed by 
Jack and Hope [ss'] which exploits memory functions and 
virtual density matrices. 

In the quantum state diffusion (QSD) method [25j . 
to obtain the operator giving the stochastic evolution 
of state vectors, one needs to perform a memory ker- 
nel integration combined with a functional derivative of 
the state vector with respect to the noise. In the NMQJ 
method the corresponding step goes in a fundamentally 
different way since the simulation produces its own non- 
Markovian quantum jump operator. This acts by trans- 
ferring the ensemble members between the existing states 
in a stochastic way [c.f. Eq. It is also worth men- 

tioning that the QSD method by definition has continu- 
ous stochastic evolution of state vectors. This means that 
in the QSD simulation iVoff ~ N. For NMQJ method, 
when the complexity of the system to be treated in- 
creases, also TVeff increases. In the ultimate limit when 
the number of different state vectors is very large, or 
even approaches infinity, then there does not exist the 
optimization scheme for NMQJ method based on Ncs- 
In this case, the simulations also become more tedious 
due to the increasing number of state vectors which need 
to evolved simultaneously. 

In general the derivation of local-in-time master equa- 



tion for driven systems is a very challenging problem in 
the theory of non-Markovian open quantum systems. We 
believe that the main difficulties here are the condition of 
very strong driving affecting the system dynamics in the 
short non-Markovian time-scale, and the case of a very 
strong coupling between the system and the reservoir. In 
the latter case, the existence of a time-local generator of 
the reduced system dynamics is not in general guaranteed 
(see section 9.2.1 of Ref. [ij). Hence, it is worth keeping 
in mind that the applicability of our method depends on 
this issue, since our starting point is the local-in-time 
master equation pT|) . 



VI. CONCLUSIONS 

We have shown that, starting from a general local- 
in-time master equation, it is possible to describe the 
dynamics of a non-Markovian open system with an en- 
semble of stochastic pure state evolutions with quantum 
jumps. The developed non-Markovian quantum jump 
method (NMQJ) demonstrates that it is indeed possible 
to unravel non-Markovian master equations with quan- 
tum jumps without making any auxiliary extensions to 
the Hilbert space of the system as done in the jump de- 
scriptions developed earlier [l^ [l^, HH [H, [2^ . Our ap- 
proach allows a rather simple and insightful description 
of non-Markovian dynamics. Even though the method 
allows to optimize the simulations in terms of using the 
effective ensemble size Nee, this number increases with 
complexity of the system under study. This sets the limit 
for the performance of the method since A'cfi state vectors 
need to be evolved simultaneously. 

The NMQJ method developed here gen eralizes a 
widely used Markovian MCWF method |l3| into the 
non-Markovian regime. Due to the existence of the 
negative decay rates for non-Markovian systems, the 
MCWF method leads to negative quantum jump prob- 
abilities. We have discovered the corresponding jump 
process which has positive probability. Due to the mem- 
ory of the system, this non-Markovian quantum jump 
essentially acts as a reverse jump, and allows the system 
to recover the information lost earlier. The consequence 
is that in the ensemble of pure states forming the den- 
sity matrix, the seemingly lost superpositions can be re- 
stored. During the time evolution, jump - reverse-jump 
cycles can occur in the ensemble members: the first jump 
during the positive decay destroys quantum superposi- 
tion while the second jump in the negative decay region 
restores them. 

Our results shed new light on the non-Markovian dy- 
namics in several ways. Breaking the density matrix 
evolution into an ensemble of state vectors with quan- 
tum jumps allows to understand how the density matrix 
carries the information on the earlier state of the sys- 
tem, and how the memory affects the system dynamics. 
This helps to clarify how local-in-time master equations 
are able to describe non-Markovian dynamics. Quantum 
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mechanics reveals often counterintuitive features. Here, 
the rate of the process appearing in the non-Markovian 
region is directly proportional to the target state of the 
process. This is opposite to the classical view where typ- 
ically the rate of a given process is given by the source 
state. Our analysis reveals in detail this counterintuitive 
feature of non-Markovian dynamics which is also present 
in the unravelled master equation. 

It has been shown earlier that Markovian open system 
dynamics with MCWF trajectories can be formally de- 
scribed as a piecewise deterministic stochastic process of 
general probability theory (36j . Consequently, we can ask 
what is the corresponding formal stochastic process for 
the NMQ J state evolutions [s^l • This holds a promise to 
exploit new stochastic process which may allow the in- 
gredients and insight by our NMQJ method to be taken 
outside the field of open systems to a more general level. 
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APPENDIX A 

In this Appendix A we show the details of the proof of 
the match between the master equation and the NMQJ 
method. 

Averaging the evolution of the ensemble 



over time step St gives 



a{t + St) = J2 



Na,{t) 

N 



j+ j- ,a' 

c,^{t)\i,^{t)){^^{t)\cUt) 



\Mt + St)){^^{t + St)\ 



J+ 



\\c,At)\Mt)W 



where we have weighted, as usual, the deterministic evo- 
lution with the no-jump probability and the jump paths 
with the corresponding jump probabilities. Above, we 
have the following quantities: P^^ {t) is the jump proba- 
bility of the state \i/ja{t)) for positive channel j+ 

PiHt) = A,4t)St{4>a.{t)\cUt)C,^{t)\^p^{t)), (A3) 



Pi->a'i^) is the reverse jump probability of state li^ait)) 
via the negative channel j_ to the state \ipa'{t)) 



Pi 



Na'{t) 



A,_it)\St 



X {^c.,{t)\Cl{t)C,_{tMo.'{t)). (A4) 



\\\Mt + st)W 

E (*)^^-«' (t)) (V-a (t) I Di-_l^, (t) 



(A2) 



r 



The reverse jump operator from the state \ipa) = 
Cj-\ipa'} /\\Cj-\^a'\\ via channel j_ to the state \ipa') 
is 

_i^^i-_,„.(i) = !V'a'(t))(^.(i)|- . . (A5) 

The deterministic evolution in Eq. (|A2p is given by 



\Mt + st)) 



iHsSt 



X |?/'a(i)), 



(A6) 



which gives for \(t)a{t + St)){(f)a{t + St)\, in first order in 
St, 



\Mt + St)){Mt + St)\ = |V'o(t))(V'a(i)| -*<5i[ffs,|V^a(t))(^a(t)|] - 2-E^J-W{^iW^^W'l^"W)(^"Wl}A^) 

J 
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It is easy to see from here, that this term gives the com- 
mutator and the anticommutator parts of the master 
equation. 

In Eq. (jA2p . the jump probabihties for the positive 
channel, appearing in the numerator and the denomina- 
tor in the no-jump path, cancel each other when doing 
the series expansion in 6t and keeping the terms to first 



order. The jump part to positive channels gives the pos- 
itive channel "sandwich term" of the master equation in 
the usual way. 

We are left with the "sandwich" term for the negative 
channels. Inserting Eqs. (|X3)) - (|T7)) into Eq. and 
comparing to Eq. ([20]) we have to show that 



^ ^\A,_{t)\stc,_mc.{t)){Mt)\cl{t) = 



N 



a a' ,j 

N, 



(A8) 



We have written here explicitly the (5-functional which 
gives the condition for the reverse jump: one can go 
via channel j- from \ipa) to \ipa') on the condition that 

IV'aW) =C',-WIV'a'(i))/l|C,_WIV'a'(i))l|- 1^ Eq. jM]), 

the last two lines cancel each other. This happens be- 
cause the (5-functional takes care of the a summation in 
the last line and the summation over a and a' are equiv- 
alent procedures making the two terms equal with oppo- 
site signs. The first and second line in Eq. (|A8p are equal. 
In the second line the (5-functional with summation over 
a means replacing \ipa) with Cj^^lipa') /\\Cj_\'ipa')\\ giv- 
ing the sandwich term of the master equation in the first 
line. Thus we have proven the equivalence between the 
master equation and the algorithm. 

The proof can be summarized in the following way: 
the deterministic part gives the commutator and anti- 
commutator parts of the master equation, the positive 
channels go in the usual way: the jump part giving the 
corresponding sandwich term of the master equation. For 
negative channels the change in the norm and jumps can- 
cel and the jump probability of negative channels times 
the deterministic evolution gives the sandwich terms. 



define short-hand notation 



APPENDIX B 



This appendix B gives the formal analytical solutions 
for the three-level systems considered in Sec. IIV CI For 
simplicity, we neglect the Lamb-shift term. First, let us 



D,{t) = / dsA,(s), 



L,{t) = / dsA,(s). 



(Bl) 
(B2) 



The direct formal solutions can be expressed by using 
these parameters, decay rates Ai{t), and initial condi- 
tions Pij(0) only. 



Two-level atom: detuned Jaynes-Cummings model 



Master equation: 



1 



I 

-iA(<){p(0,(T+a_}. (B3) 



Jump operator: 



Ci =a_ = |6)(a|. 



(B4) 



Populations: 

Paait) = e-^i(*Vaa(0), (B5) 
Pbbit) = {l-e-^l(*)}paa(0)-f Pb6(0). (B6) 

Coherences: 

Patit) = e-^^(*)/Vah(0). (B7) 
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Three-level atom: A-system 



Coherences: 



Master equation: 



+ Ai(t) 

+ A2(f) 



\b){a\p{t)\a){h\-\{p{t),\a){a\} 
\c){a\p(t)\a){c\-\{p{t),\a){a\} 



(B8) 



Jump operators: 

Ci = |6)(a|, (B9) 
C2 = \c){a\. (BIO) 

Populations: 

p,,(t)=e-[^^W+^^(*)Vaa(0), (Bll) 
Pbbit) = /*dsAi(s)e-[^^(^)+^=(^)lpaa(0) 

Jo 

+ PbbiO), (B12) 

Pec(i)= / dsA2(s)e-[^^W+^=Wlpaa(0) 
JO 

+ Pcc(0). (B13) 

Coherences: 

= e-['^i(*)+*^=(*)+^i(*)/2+C2(t)/2]^^^(0)^ (B14) 
= e-[^^i(*)+^i2(*)+^i(*)/2+O2(t)/2]^^^(0)^ (B15) 
Pbc{t) = Pbc{0). (B16) 

Three-level atom: F-system 

Master equation: 

m =]x,ma){a\,p{t)] + -^x2mb){b\,p{t)] 



+ Ai(t) 

+ A2(i) 



|c)(a|p(t)|a)(c|--{p(t),|a)(a|} 



c){b\p{t)\b){c\--{p{t),\b){b\} 



Jump operators: 

Ci = \c){a\, 
C2 = \c){b\. 

Populations: 

Paait)=e-'''^'^Paa{0), 

Pbbit) = e-^=Wp66(0), 



Pccit) 



+ Pcc(O). 



Paa(O) + 



1 - e-^=W 



(B17) 

(B18) 
(B19) 

(B20) 
(B21) 

PbbiO) 
(B22) 



Pabit) 
Pacit) 
Pbcit) 



^ -[iLi(t)+iL2(t)+Di(t)/2+D2it)/2] 



^ g-[jLi(t)+_Di(t)/2] 
-[iL2{t)+D2{t)/2 



PaciO), 

PbM- 



PabiO), (B23) 
(B24) 
(B25) 



Three-level atom: Ladder-system 



Master equation: 



Pit) =-^X^it)[\a){a\,pit)] + h2it)[\b){b\,pit)] 



+ Ai(i) 

+ A2(t) 

Jump operators: 



\b){a\p{t)\a){b\-Up{t),\a){a\} 



\c){b\pit)\b){c\--{p{t),\b){h\} 



(B26) 



Ci = \b){a\, 
C2 = \c}{b\. 



(B27) 
(B28) 



Populations: 



Paait) = e-^i(*Vaa(0), (B29) 

Jo 



-D2{t), 



e -^-^>PbhiO), 



(B30) 



Pccit) 



1 - e 



/ dsAi(s)e-^^W+^^Wlp„„(0) 
Jo J 



1 - e 



P6fc(0)+Pcc(0). (B31) 



Coherences: 



= e-I'^iW-'^^W-^iW/2-^^(*)/21p„b(0), (B32) 
Pac(i) = e-I'^^W+^^W/^lpac(0), (B33) 
p,,(f) = e-['^=W+^^W/2V6c(0). (B34) 
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